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Abstract 

Chaos transition, as an important topic, has become an active research subject 
in non-hnear science. By considering a Dicke Hamihonian coupled to a bath of 
harmonic osciUator, we have been able to introduce a logistic map with quantum 
corrections. Some basic dynamical properties, such as Lyapunov exponents and 
bifurcation diagram of the model are studied, we show that in this model, the 
transition from integrable motion to periodic, chaotic and hyperchaotic as the 
control parameter r is increased. 
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1 Introduction 



Investigation in the area of quantum dots couple with another quantum dot or one 
external field has received a great deal of attention in recent years. Deterministic equa- 
tions have been used nowadays in many applications [1] . For instance, a chaotic system 
is one which displays complex and unpredictable behavior. These systems with some 
of properties such as sensitivity to initial condition [2j, system parameter, ergodicity 
[3] and mixing |^]. One of the popular fields of quantum chaos addresses quantum 
maps. Classical maps exhibit the properties of chaotic dynamical systems. Various 
chaotic maps use in Cryptography, Watermarking, Random number generators and 
lots of other fields. On the other hand, a big number of researches have been dedicated 
so far to quantum maps as paradigms of quantum chaos [5l|6]. M.Goggin et al. in their 
work, the quantum correction were made by introducing a periodic noise to a model of 
oscillatory system coupled to a bath of oscillators. They derive a new logistic map with 
quantum computing that illustrate a period-doubling route to chaos The quan- 
tization schemes for maps are different and many classically chaotic maps have been 
quantized including the standard map [H [9], logistic map [7], baker map on the torus 
[TOl im [T2] and on sphere [H]. In present paper we introduce generahzed Dick model 
has been considered to present a new quantum map which demonstrates positive lya- 
punov exponents. The comparison between analytically studied and simulation results 
indicate that the propose methods provide good estimates. The paper is organized as 
follows. In section 2, introduces the model. In section 3, we calculate the lyapunov 
exponents and KS-entropy. Finally, in section 4, we conclude the paper. Respective 
appendix can be found at the end of paper. 



H = a)a + ujaJz + -^{a^ + a){J^ + J+) + — =1/( J_, J+)S„5(t - nT), (1) 



this hamiltonian is the generalized Dick Hamiltonian. In fact we add a delta function 
to Dick Hamiltonian. Therefore, a and a) bosonic annihilation and creation operators, 
energy separation of N two- level atoms equal to tuJjA , h = 1 , field of frequency Uf, 
cj^ = cI'a/'^/ > is given in units of the frequency of the field, 7 = 7/0;/ is the coupling 



2 The model 



We start with the hamiltonian: 
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parameter. the atomic relative population operator, and J± the atomic transition 
operators. The commutation relations for the fermionic operator J are: 

[J+, J_] = 2 J, , [J„ J±] = ± J± , [a, at] = 1. 

Heisenberg equation of motion reads: 



by considering (-i^[J+, y( J+, J_)] = /( J_, J+)), 
we have : 

J+ = iujaJ, - i^{a^ + a)(2J,) + /(J+, J_)S„(5(t - nT), {Al) 



(2) 



d = -m-i^(J_ + J+) 



(^2) 



In order to obtain an equation for J+(i), we substitute the solution of (A2) into 
(Al): 

J+ = iujAJ+ - i^at(0)(exp(it))(2J^) - i^a(0)(exp(-it))(2J^) 
V v A^ 



+/(J+,J_)E„5(i-nr), (3) 
we integrate of Eq. (3) from t= nT to t= (n+l)T, so as to obtain 



J+((n + l)r) - J+{nT) ex.-p{iujAT) = 

(n+i)r 



r (n+l)T 

-i^=aUQ) / (2J^ exp(i(l - ujA)t)dt 

^ JnT 



z^a(O) / (2J2)exp(-«(l + a;A)t)c^t 

A" JnT 



/(J+(nT),J_(nT))exp(^c^^T). (4) 



We assume that < a(0) >= , < 0^(0) >= 0, the expectation value of J+ can be 
expressed in: 



< J+((n + 1)T) >= (< J+{nT) > + < f{J+{nT), J^nT)) >) eicp{iu;AT). (5) 

For comparison our quantum results to the familiar logistic map we choose the 
"force": 



f{J+{n)^ J -in)) = -J+{n) + f (■^+(n) " J~{n)J+{n)) QW{-i^AT), 



(6) 



3 



where r is an adjustable parameter. Eq. (5) then becomes : 

< J+(n+l) >= r (< J+(„) > — < J-{n)J+(n) >) , (7) 

where x =< J+ >=< J_ > , then: 

Xn+l l^i^Xn Xn )• 

Which is the classical logistic map . 

The logistic map is one of the most studied discrete chaotic maps. It was first 
proposed as pseudo random number generator by Von Neumann in 1947 partly because 
it had a known algebraic distribution and mentioned later, in 1969, by Knuth[T4| [15]. 
The logistic map is given by: 

Xn+l = rXn{l — Xn) (8) 

where Xn G (0, 1) and r are the system variable and parameter, respectively, and n is 
the number of iterations. Thus, given an initial value xo and a parameter r. In this 
paper, we refer to xq and r as the initial state of the logistic map. 

Xn+i = rXn{l — Xn), foT Xn € (0, 1), and r G (3.99996,4] 

and its behavior is full chaotic(Fig 1). The choice of r in the equation above guarantees 
the existence of a chaotic orbit that can be shadowed by only one map as stated in |16] . 
The lyapunov exponent h is taken as a measure of the "chaoticness" of N quantum dots 
coupled with a single bosonic mode. The dynamical systems, even the discrete-time and 
one-variable ones have different types of behaviour. The system can be in a fixed point 
and nothing changes, the trajectory of the system may also be a cycle with a certain 
period. Fixed point and periodic orbits may be stable or unstable. We are usually in- 
terested in an invariant measure , i.e. a probability measure that does not change under 
the dynamics. The probability measure p{x) = — , ^ on [0; 1] is a Sinai-RuUe-Bowen 

(SRB) measure which is an invariant measure which describes statistically the station- 
ary states of the system and absolutely continues with respect to Lebesgue measure. 
Entropy sometimes is called uncertainty and sometimes information. In order to study 
the stability, entropy can be used as an acceptable parameter. Also, by considering SRB 
measure, it is possible to continue this study with KS-entropy, the well-known measure 
for chaos in dynamical system. This section introduces the KS-entropy for hierarchy of 
pair-coupled maps with dynamical parameter. We try to calculate Lyapunov exponent 
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as another tool to study the stabihty. 



h= An \F'{x)\p{x)dx = 2 A'y^lLJi^idx, (9) 

Jo Jo 'Ky/x{l — x) 

with change of variable x = sin^(6') , dx = 2sin(^) cos{9)d9, and using Eq. (9) becomes 

r^/Hn\r{l-2sm\e))\^ 



Jo 7rsin(&j cos(&) 



(10) 



KS-entropy, may also can be written as: 



/i = lnr-ln2, (11) 

for r = 4 — > /i = In 2 > 0, it is positive. 

For appearance the effect of the quantum correlations by using J+ =< J+ > +SJ+ 
and J_ =< J_ > +SJ-, one gets 

< J+(n+i) >= r (< J+(„) > - < J_(„) >< J+(„) >) - r < 5J_5J+ > . (12) 

First, the time derivative of(5J+5J_) is 

^{dJ+dJ_) = dJ+dJ_ + dJ+dJ_. (13) 

Inserting these two expressions J+ =< J+ > +5J+ and J_ =< J_ > +(5J_ into Eq. (6) 
we end up with 

/( J+, J_) = - < J+ > -5J+ + exp(— i-u;yiT)r(< J+ > < J_ >< J+ > - 

6J_6J+- < J_ > 5J+ - 5J_ < J+ >). (14) 

Inserting J+ = and — SJ^ , into Eq. (3) we obtain: 

SJ+ = iwA{< J+ > +5J+) - i^at(0)(exp(it))[5J+5J_ - 5J_5J+] - i^a(O) 



(exp(-ii))[5J+5J_ - 5J-5J+] + [- < J+ > -5J+ + exp(-iwAr)r(< J+ > 

- < J_ >< J+ > -SJ^SJ+- < J_ > 5J+ - 5J_ < J+ >)]^J{t - nT). (15) 
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Due to 



5J_ = 5J+^ = -iwA{< J- > +SJ-) + i-J=a{0){exp{-it))[6J_6J+ - 6J+SJ_] + i--^a^(O) 

(ex.p(it))[5J^5J+ - 5J+5J^] + [- < J_ > -5J_ + exp(iwAT)r(< J_ > +5J_ - 

< J+ >< J_ > -5J+SJ^- <J+> 5J_ - SJ+ < J_ >)]J:j{t - nT),{16) 

Finally, inserting Eq. (15) and Eq. (16) into Eq. (13), we end up with 



iwA{< J+ > +SJ+) - i^a^{0){exp(U))[SJ+SJ_ - 5J_5J+] - i^a(O) 

y N y N 



(exp(— it))[(5J+(5J_ - 5J_(5J+] + [- < J+ > -SJ+ + exp{-iwAT)r{< J+ > +6J+- < J_ > 
< J+ > -6J_6J+- < J_ > 6J+ - 6J^ < J+ >)] EJit - nT)] 5J_ + 6J+ [-iwA{< J- > +SJ-) 
+i^a{0){exp{-it))[5J-SJ+ - 5J+5J-] + i^a\0){exp{it))[5J-5J+ - 5J+5JJ\ + [- < J_ > 



-5 J- + exY>{iwAT)r{< J_ > +5J_- < J+ >< J_ > -SJ+SJ^- < J+ > SJ- - 5J+ < J_ >)] 

EJ{t - nT)] .(17) 



By integrating Eq. (17), from (nT) to ((n+l)T), and take the expectation value, by 
taking into account < 5J-{nT) >=< 5J+{nT) >= , < a^{0) >=< a(0) >= we 
obtain: 



< SJ+{{n + l)T)SJ-{{n + 1)T) >= - < SJ+{nT)5J-{nT) > +rexp{-iwAT) < SJ+{nT)5J-{nT) > 
-rexp{-iwAT) < J_ >< 6J+{nT)6J_{nT) > -rexp(-iw^T) < 6J_{nT)6J_{nT) >< J+ > +r 
exp(iwyiT) < 5J^{nT)5J_{nT) > —rexp{iwAT) < 5J^{nT)5J_[nT) >< J_|_ > — rexp(i'u;^T) 

< 5J+{nT)5J+{nT) >< J_ > . (18) 



The calculation of < > goes as follows: 



^{5J+5J+) = 5J+5J+ + 5J+6J+. 



(19) 



appearance the effect of the quantum correlations by using J+ =< J+ > +(5J+ and 
J_ =< J_ > +(5J_ Eq. (6), one gets: 

/( J+, J_) = - < J+ > -5J+ + exp(-iM;^T)r(< J+ > +5J+- < J_ >< J+ > - 

5J-6J+- < J_ > 5J+ - 5J_ < J+ >), (20) 

Inserting = 6J^ , into Eq. (3) we obtain: 



SJ+ = iwA(< J+ > +SJ+) - i^a\0){eicp{it))[SJ+SJ^ - SJ_SJ+] - i^a(O) 



{exp{-it))[6J+SJ_ - 5J_5J+] + [- < J+ > -5J+ + exp(-i;w^T)r(< J+ > 

- < J_ >< J+ > -6J_6J+- < J_ > 5J+ - 5J_ < J+ >)]S„5(t - nT), (21) 

Finally, inserting Eq. (21) into Eq. (19), we end up with: 

^{SJ+SJ+) = [iwA{< J+ > +5J+) - i^a^(0)(exp(it))[5J+5J_ - 5J_5J+] - i^a(O) 
at y/N vN 



(exp(-it))[(5J+5J_ - 5J_5J+] + [- < J+ > + cxp(-^^^;^T)r(< J+ > +5J+- < J_ > 

< J+ > -6J-6J+- < J_ > - 5J_ < J+ >)]E„5(t - nT)]5J+ + SJ+[iwA{< J+ > + 
SJ+) - i^aU0)(exp(it))[SJ+SJ- - SJ-SJ+] - i^a(O) 



{exp{-it))[SJ+SJ- - SJ-SJ+\ + [- < J+ > -6J+ + exp{-iwAT)r{< J+ > +6J+- < J_ > 

< J+ > -SJ_SJ+- < J_ > SJ+ - 6J_ < J+ >)]E„5(t - nT)]. (22) 

By integrating Eq. (22), from (nT) to ((n+l)T), and take the expectation value, 
by taking into account < 5J-{nT) >=< 5J+{nT) >= , < a^(0) >=< a(0) >= we 
obtain: 

< SJ+{{n+l)T)SJ+{{n+l)T) > exp{-2iiUA{n + 1)T)- < SJ+{nT)SJ+{nT) > exp{-2iiUAnT) 
exp{-2iujAnT){- < 5J+{nT)5J+{nT) > +exp{-iujAT)r{< 5J+{nT)5J+{nT) > - < J^{nT) > 
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< dJ+{nT)SJ+{nT) > - < J+{nT) >< 5J_{nT)5J+{nT) >)) + exp(2ia;Anr)(- < 5J+{nT) 
5J+{nT) > + exp{-iUAT)ri< 5J+{nT)5J+{nT) > - < J_{nT) >< SJ+{nT)6J+{nT) > - 



By taking into account < SJ+{nT)5J_{nT) >=¥„ , < >—Zn , < J+{nT) >—X, 

and inserting these values into Eqs. (12), (18) and (23), we obtain: 



Y(„+i) = -Fn + rexp(-iwAT)([(l - X^^)* 

+ exp(2iw^r) - exp{2iwAT))Y^] - - exp(2iwAr)Z„X„*) {B2) 

^ Z{n+i) = -Zn exp{2iwAT) + r exp{iwAT){2Zn - 2X„*Z„ - - F„X„) {B3). 

These equations give us lowest-order quantum corrections.We consider (3 — iwaT. 

3 The onset of quantum chaos 

The research on the chaos of Eq. (B) can begin with the analysis on the stability of the 
fixed point. Assume (x, y, z) is the fixed point of Eq. (B), then (x, y, z) is the solution 
of the equations below: 

X — r{x — x^) — ry (CI) 



< SJ+{nT)6J_{nT) >< J+{nT) >)). 



(23) 



{Bl) 



y = -y + rexp{-/3){[{l - X 

+ exp(2/3) — X exp{2(3))y] — zx — exp{2(3)zx) 



(C2) 



z = —z exp{2/3) + r exp{/3){2z — 2xz — xy — yx) (C3). 
By considering exp(/3) = a and solution Eq. (C) we obtain: 
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X = 0, 
(r - l)/r, 

— (— *r + 2*a^*r^ + 2*a^ — 6*a^*r + 2*a*r^ + 2*a — r)/(a^*r — 4*a^*r^ + 
6*a^*r — 4*a*r^ + r), 

y = o, 

0, 

(— a^*r — 2* a''* + 6* a'^*r^ + 2* a^*r + 2* a''' + 4* a^*r'^ — 8* a®* r'^ — 20 *a^*r — 
4 * a'^ - 14 * * + 42 * * + 14 * * r + 14 * + 8 * 0,^ * r"^ - 16 * a"^ * - 54 * a'^ * r - 
8 * — 14 * * + 42 * * + 14 * * r + 14 * + 4 * * r'^ — 8 * * — 20 * * r — 
4*a^ — 2*a*r^ + 6*a*r^ + 2*a*r + 2*a — r)/(a^*r^ — 8* a'^*r^ + 16* a®* r"^ + 12* a^* 
r^ — 56* a^*r^ + 32 *a^*r^ + 38* a^*r^ — 56*a^ *r^ + 16* a^*r'^ + 12 *a^*r^ — 8*a*r^+r^), 

^ = 0, 
0, 

(2 * * + 4 * * r"^ — 8 * * — 4 * * — 8 * * r — 16 * a * r''^ + 48 * a * + 16 * a '' * 
r + 8*a'^ + 12*a^*r^-2A*a^*r^ + 4*a^*r^-88*a^*r-16*a^' -32*a^*r^ + 92*a^* 
r2 + 32*a^*r + 48*a^ + 12*a^*r^-24*a^*r3 + 4*a^*r2-88*a^*r-16*a^-16*a3* 
+ 48 * * + 16 * *r + 8*a^ + 4* a^*r^ — 8* a^*r^ — 4* a^*r^ — 8* a^*r + 2*a* 
r^)/(a^O * - 8 * * + 16 * * + 13 * * - 64 * * + 48 * * + 50 * * - 
112 * *r^ + 48* *r^ + 50 * *r^ - 64 *a^*r^ + 16* a^*r^ + 13 *a^*r^- 8* a *r^ + r^). 

In the thermodynamic hmit /3 — > or a = 1, we have following fixed points: 

(xo,yo,^o) = (0,0,0) 
(xi,yi,zi) = ((r-l)/r,0,0) 

(2^2, y2, -^2) = (-(r^ - 3 * r + 2)/(4 * r - 3 * r^) , (2 * r'' - 10 * + 18 * - 14 * r + 4)/(9 * 
- 24 * + 16 * r^), (r^ - 6 * + 13 * - 12 * r + 4) /(9 * - 24 * + 16 * r^)) 



We substitute X(n+i) = X, X„ = x, y(„+i) = F, = y, Z(^n+i) = -Z' and ^(n+i) = z 
into Eq. (B) obtain: 
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X = r{x — x"^) — ry 



Y 



-y + r cxp(— /?)([( 1 — X 



+ exp(2/3) — X exp{2/3))y] — zx — exp{2/3)zx) 



m 



Z ^ -z exp{2f3) + r exp{P){2z - 2xz - xy - yx) (D3). 

The corresponding bifurcation diagram of state x with respect to r is giving in Fig. 2. 

In order to calculate the lyapunov exponents at x — Xq, y — yo and z — Zq, we need 
to calculate the characteristic roots of the matrix: 



- dX 


ax 


dx - 




dx 


dy 


dz 




9Y_ 


dY_ 


dY 




dx 


dy 


dz 




dZ 


dZ 


dZ 




. dx 


dy 


dz . 





r(l - 2xo) 
l(yo + zo)(-l-a') 
2ar{-Zo - yo) 



—r 

-l + ^(l-xo)(l + a2) 
—2raxo 



Substituting xq — 0, yo — and = , into A: 



A 



-l + i(l + a2) 




and obtain three eigenvalues as below: 






-a^ + 2ra 





X-1 - a') 
-0? + 2ra(l — Xo) 
(24) 



(25) 



r — A 





-r 
-l + ^(l + a2)-A 

-a^ + 2ar - A 



Ai = r, A2 = -1 + -(1 + a^), A3 = -a^ + 2ar 
a 



= (r-A)(-l+-(l+a2)-A)(-aV2ar-A) = 
a 



(26) 
(27) 



for a = 1, Ai = r, A2 = -1 + 2r, A3 = -1 + 2r. 
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As increased r, the stability of system (D) is summarized as follows: 

(1) r = 0, Ai = 0, A2 = A3 < 0, system (D) is periodic. 

(2) < r < 0.5, Ai > 0, A2 = A3 < 0, the stability of system (D) is chaotic. 

(3) r > 0.5, Ai > 0, A2 = A3 > 0, the stability of system (_D) is hyperchaotic [T7] . 



Obviously, the coupled system is ergodic at (x = xo,y = yo, z = zo) according to 
Pessin theorem [181 119]. Due to ergodcity of one-dimensional map, Xn+i = — 
we have Ai, A2, A3, and the KS-entropy is equal to sum of positive Lyapunov exponents. 
Now comparing the KS-entropy with sum of Lyapunov exponents, one can show that: 



h 



KS 



Ai + A2 + A3 = 5r - 2. 



(2^ 



For r = 0.4 

hKS = 0, 

therefore the system is not chaotic. The measurable dynamical system is chaotic for 
hxs > and predictive for hxs = 0. For r > 0.4, positive KS-entropy occurs. For each 
r > 1, set A2 = A3 > Ai. By considering 0.5 < r < 1 ^ Ai > A2 = A3 > 0. 

The corresponding bifurcation diagram of state x with respect to r is giving in Fig. 
3. Substituting x = — , y = and z = we into A: 



— r + 2 — r 
A= -1 + 1 + a 
-2a(r-l) 

and obtain three eigenvalues as below: 





r-l 
a 

-a^ + 2a 



(29) 



-r + 2 - A 





-r 

-1 + i + a-A ^ 

a a 

-2a{r - 1) -a^ + 2a - A 



-r+2-A)(-l + -+a-A)(-a2+2a-A) = 
a 



Ai = -r + 2, A2 = -IH h a, A3 = -a^ + 2a 



a 
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(30) 



for a = 1, Ai = 2 — r, A2 = 1, A3 = 1. 

As increased r, the stability of system (D) is summarized as follows: 

(1) r < 2, Ai > 0, A2 = A3 = 1 > 0, system (D) is hyperchaotic. 

(2) r = 2, Ai = 0, A2 = A3 = 1 > 0, the stability of system (D) is hyperchaotic. 

(3) r > 2, Ai < 0, A2 = A3 = 1 > 0, the stabihty of system (D) is hyperchaotic. 

Now comparing the KS-entropy with sum of Lyapunov exponents, one can show that: 

/iKS = Ai + A2 + A3 = 4 - r. (31) 

For r — 4: 
hxs = 0, 

therefore the system is not chaotic. The measurable dynamical system is chaotic for 
hxs > and predictive for Hks = 0. For r < 4, positive KS-entropy occurs. For each 
r > 1, set A2 = A3 > Ai. By considering r < 1 =^ Ai > A2 = A3 > 0. 
The corresponding bifurcation diagram of state x with respect to /3 is giving in Fig. 5. 
In accordance with this figure, the system is chaotic for positive (3. 

4 Conclusion 

In this article, a new scheme for studying the effects of quantum correlations on a Dicke 
Hamiltonian coupled to a bath of harmonic oscillator has been presented. By calculating 
Some basic dynamical properties, such as Lyapunov exponents and bifurcation diagram 
of the model, the proposed scheme is proved to be the transition from integrable motion 
to periodic, chaotic and hyperchaotic as the control parameter r is increased. 

Further, we do hope that our obtained results through this paper will pave the way 
for further studies on quantum chaos. 

5 Appendix 

Here we derive the Lyapunov exponents of chaotic maps by using fixed point (x2, y2, ^2)- 
For a = 1, we have following fixed points: 
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(x2,y2,Z2) = (-(r^ - 3*r + 2)/(4*r-3* r^), (2 * - 10 * + 18 * - 14 * 
r + 4)/(9*r^-24*r^ + 16*r2), (r^-6*r^ + 13*r2-12*r + 4)/(9*r^-24*r^ + 16*r2)) 

for r = 4, 

(2^2,^2,-^2) = (0.19,0.105,0.352). 

The corresponding bifurcation diagram of state x with respect to r is giving in Fig. 
4. Substituting x = 0.19, y = 0.105 and z = 0.352 we into A: 



2.48 -4 

-3.656 5.48 -1.52 

-3.656 -1.52 5.48 

and obtain three eigenvalues as below: 



(32) 



2.48 - A 
-3.656 
-3.656 



-4 
5.48 - A 
-1.52 





-1.52 
5.48 - A 



-A^ + 13.44A2 - 40.276A - 33.618) = 



Ai=24919/(3750*(-2138981/125000+(10298398109(l/2)*16875000000(l/2)*i)/16875000000)(l/3))+ 

(-2138981/125000 + (10298398109^1/2) * 16875000000^1/2) * i)/16875000000)(l/3) + 

112/25 

A2=112/25-(-2138981/125000+(10298398109(l/2)*16875000000(l/2)*i)/16875000000)(l/3)/2- 
24919/(7500*(-2138981/125000+(10298398109(l/2)*16875000000(l/2)*i)/16875000000)(l/3)) + 
(3(1/2) * (24919/(3750 * (-2138981/125000 + (10298398109^1/2) * 16875000000^1/2) * 
i)/16875000000)(l/3)) - (-2138981/125000 + (10298398109^1/2) * 16875000000^1/2) * 
i)/16875000000)(l/3)) *i)/2 

A3=112/25-(-2138981/125000+(10298398109(l/2)*16875000000(l/2)*z)/16875000000)(l/3)/2- 
24919/(7500*(-2138981/125000+(10298398109(l/2)*16875000000(l/2)*i)/16875000000)(l/3))- 
(3(1/2) * (24919/(3750 * (-2138981/125000 + (10298398109(1/2) * 16875000000(1/2) * 
i)/16875000000)(l/3)) - (-2138981/125000 + (10298398109(1/2) * 16875000000(1/2) * 
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i)/16875000000)(l/3)) *i)/2 
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Figure 1: Bifurcation diagram for the logistic map, Xn+i = tx„(1 — x„) when the 
parameter r is varied from 1 to 4. 




0.5 1 1.5 2 2.5 3 3.5 4 



Figure 2: Bifurcation diagram for increasing r, with a = 1, xq = 0, yo = and = 0. 
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Figure 3: Bifurcation diagram for increasing r, with a = 1, xq = yo = and zq = 0. 
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Figure 4: Bifurcation diagram for r = 4, with a = 1, xq = 0.19, yo = 0.105 and 
2o = 0.352. 
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Figure 5: Bifurcation diagram for increasing (3, with r = 4, xq = 0.2, yo = and 
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